function [err] = SC(DFILE, sf, plot_flag)
% Single channel data analysis.
% source taken from VC_IV 9 Dec 1999. P. Manis
% Calculate IVs based on data passed through DFILE, CONTROL(sf), current and voltage
% plot_flag : 1 is to display, 2 is to print, and 3 is to print and then close
% Data is assumed to have been read and corrected when it is passed to this routine
% Step times for for the data are derived from contents of the CONTROL matrix
% rather than the DFILE values, so we can correct them if necessary.

global VOLTAGE, CURRENT
global CONTROL

err = 0; % save output in case we exit early...
protocol=lower(CONTROL(sf).protocol);
QueMessage(sprintf('Processing %s, %s for IV', DFILE.filename, protocol));
rate = DFILE.rate*DFILE.nr_channel/1000;
[records,pts]=size(current);
time=make_time(DFILE);
tmax=max(time);


% mean times for standard IV are computed as percentages from the #2 window
ts1=number_arg(CONTROL(sf).durho);
ts2=number_arg(CONTROL(sf).durs1);
ts3=number_arg(CONTROL(sf).durs2);
ts4=number_arg(CONTROL(sf).durs3);

% analysis mode definitions.
HIST=1;
PofV=2;

analysis=0; % none set.

protocol=lower(deblank(strjust(protocol, 'left'))); % remove leading and trailing spaces from the string.
% set up parameters for analysis
switch(protocol)
case 'pativ' % standard patch IV from -60 mV, or with prepulse
   analysis=bitor(HIST, PofV);
   tdelA=ts1; tstepdurA=ts2; % one step, right after holding.

case {'patqp2' 'qinactp2'} % patch IV with prepulse 
   analysis=bitor(HIST, PofV);
   tdelA=ts1+ts2; tstepdurA=ts3; % one step, right after holding.

case 'ivvar' % for variance analysis
	analysis=bitor(HIST, PofV);
	tdelA=ts1; tstepdurA=ts2;
case 'hold' % patches just held a various voltages across the traces
   analysis=bitor(HIST, PofV);
   tdelA=0; tstepdurA=ts1+ts2+ts3+ts4; % set window
otherwise
   QueMessage(sprintf('SC: processing for protocol %s not implemented', upper(protocol)));
   return;
end

% compute windows for activation analysis
if(bitget(analysis,HIST)) % windos for histogram analysis 
   hist_1=tdelA+tstepdurA*0.01;
   hist_2=tdelA+tstepdurA*0.99; % most of window - but not ends
   
   h_t1=floor(hist_1/rate);
   h_t2=floor(hist_2/rate);
end

if(bitget(analysis,PofV)) % window for open probability analysis
end

%
% filter the current traces
% According to Colquhoun and Sigworth (Pratical Analysis of Records, Chapter 19 in
% Sakmann and Neher, Single Channel Recording)
% the optimal filter is Gaussian for the first
% pass of analysis. So, we use here their Gaussian filter.
% the filter itself is in Gaussian.m in pbm_tools

fsamp = 1000/rate; % get sampling frequency
fco = 0.1;		% cutoff frequency units of sample frequency
if(fco < 1) % if fco is > 1 then this is not a filter!
   current = gaussian(current, fco); % filter all the traces...
end

% now compute various parameters.
% the relevant IVs.
% there are up to 4 possible IVs:
% volt_ssA is the activation voltage; volt_ssI is the voltage for inactivation protocols.
% curr_peak_A is peak current for activation protocols; curr_ss_A is the peak current
% curr_peak_I is peak current fir inactivation protocols; curr_ss_I is the steady state current

if(bitget(analysis, HIST)) % compute histograms for every voltage level in the set.
   QueMessage('Computing histograms');
   % first we have to find the voltage levels...
   delta=2.5; % delta for voltages, in mV.
   nbins=101; % use odd number...
   xv=[-150:delta:150]; % voltage levels to test
   mv=mean(VOLTAGE(:,h_t1:h_t2)'); % voltages in data set
   [n, mo]=hist(mv, xv); % makes a histogram of the Voltages...
   w=find(n>1); % find bins with an entry
 
	mn=min(min(CURRENT(:,h_t1:h_t2)));
	mx=max(max(CURRENT(:,h_t1:h_t2)));
	outbin=[mn:(mx-mn)/(nbins-1):mx];
   k = 0;
	sh=zeros(1,nbins);
   vh=0;
	rl.r=[];
	rl.v=[];
	for i=w
   	  [z a]=find(abs(mv-xv(i))<= delta); % find all records corresponding to that data
	  if(~isempty(z)) % we have data in the range.
	  	[h1, ix] = hist(CURRENT(a,h_t1:h_t2)',outbin);
		k=k+1;
		sh(k,:)=sum(h1'); % sum across trials.
		vh(k)=xv(i); % the nominal value is here.
	  rl(k).r=a; % records (offset...)
	  rl(k).v=xv(i); % and corresponding voltage
		end;
   end;
   nhist=k;
	% sort record list by voltage.....
	[r,u]=sort([rl.v]); % is already in order, but lets be sure..
	rvlist=[rl(u).r]; % that's it.	
	sh=[sh(u,:)];
	QueMessage('Complete computing Histograms');
end

% ctl=CONTROL; % return the modified structure with analysis results

% now create the figure

downsample=1;
fsize=7;

if(~exist('plot_flag')) % if there is no argument, just do the plot.
	plot_flag=1;
end;
if(~plot_flag)
	return;
end;
h = findobj('Tag', 'SC'); % check for pre-existing window
if(isempty(h)) % if none, make one
   h = figure('Tag', 'SC', 'Name', 'Single Channel Analysis', 'NumberTitle', 'off');
end
figure(h); % otherwise, select it
clf; % always clear the window...
set(gcf, 'Backingstore', 'off');
orient landscape

%plot currents and stack traces
% first, we reshape the data into groups, end-to-end
g=CONTROL(sf).groupavg; % this determines the size of a group
cu = CURRENT;
cu(:,(h_t1-1))=NaN;
cu=cu(rvlist, (h_t1-1):h_t2);
[Nrec, Npoints]=size(cu);
if(mod(Nrec, g) ~= 0) % whoops, we don't have a full deck to work with. Take it up to the next level
	Nadd=g*ceil(Nrec/g)-Nrec;
	ada=NaN*zeros(Nadd, Npoints); % add a few col's of NaN's (no data...)..
	cu=[cu; ada]; % we can do this to augment the array (also to add the average...)
	[Nrec, Npoints] = size(cu);
end;
cu=reshape(cu', Npoints*g, Nrec/g)'; % now we should be safe.
[NNewr, NNewp] = size(cu);
pl_min=min(min(cu))/(NNewr*1/3);
pl_max=max(max(cu))/(NNewr*1/3);
pl_y=max(abs(pl_min), abs(pl_max));
% make a matching time base 
h_sp1=subplot('Position',[0.07,0.1,0.65,0.85])
%make time array

ioff=pl_y*([1:NNewr]-NNewr/5);
for i=1:NNewr
	cu(i,:)=cu(i,:)+ioff(i);
end;
% make the time base
tbb=time((h_t1-1):h_t2)-time(h_t1); % starts a 0
tbb(1)=NaN;
tbmax=1.05*max(tbb);
tb=[];
for i=1:g
	tb= [tb tbb+(i-1)*tbmax];
end;
line(tb,cu, 'Parent', h_sp1, 'EraseMode', 'none', 'LineWidth', 0.1);
set(gca, 'FontSize', fsize);
ylabel('nA')
u=get(gca, 'YLim');

%plot voltage commands
%h_sp2=subplot('Position',[0.07,0.05,0.65,0.15]);
%line(time,voltage, 'Parent', h_sp2, 'EraseMode', 'none');
%set(gca, 'FontSize', fsize);
%ylabel('mV');
%xlabel('ms');
%set(gca, 'Xlim', v);

% plot histograms at each voltage...
dy=0.9/(nhist+1);
for i=1:nhist
	h_hsp(i)=subplot('Position',[0.75,0.95-dy*i,0.2,dy]);
	bar(outbin,sh(nhist-i+1,:));
	set(gca, 'FontSize', fsize);
	if(i < nhist)
		set(gca, 'XTickLabel', '');
	end;
	xm=get(gca, 'XLim');
	ym=get(gca, 'YLim');	
	text(xm(2),ym(2),sprintf('%6.0f mV',vh(nhist-i+1)),'Parent', h_hsp(i), 'HorizontalAlignment', 'right', 'VerticalAlignment', 'top', 'FontSize', 6);
	ylabel('N');
	if(i == nhist)
		xlabel('pA');
	end;
end;	

%plot text stuff
h_info=subplot('Position',[0.5,0.0,0.25,0.05])
r_uncomp = (1-number_arg(CONTROL(sf).comp/100))*number_arg(CONTROL(sf).access);
axis([0,1,0,1])
axis('off')
text(0,0.75,sprintf('%-12s R[%d:%d]   %-8s',DFILE.filename, DFILE.frec, DFILE.lrec, CONTROL(sf).protocol), 'Parent', h_info, 'Fontsize', 8);
text(0,0.5,sprintf('Solution:%-12s  gain:%4.1f  LPF:%4.1f kHz', CONTROL(sf).solution, DFILE.igain, DFILE.low_pass(1)), 'Parent', h_info, 'FontSize', 7);
text(0,0.1,sprintf('%s',date),'Parent', h_info, 'FontSize',6);

orient landscape

return;
if plot_flag>1
   print
end
if plot_flag==3
   close
end
